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Abstract 

A  review  of  recent  literature  on  proton  exchange  membrane  fuel  cell  modeling  is  presented.  Fuel  cell  models  can  be  categorized  as 
analytical,  semi-empirical  or  mechanistic.  Mechanistic  models  can  be  further  subcategorized  based  on  the  solution  strategy,  single-domain  or 
multi-domain.  The  multi-domain  approach  develops  and  solves  separate  equations  in  each  region  of  the  fuel  cell.  The  single-domain  approach 
consists  of  equations  governing  the  entire  domain  of  interest,  with  source  and  sink  terms  accounting  for  species  consumption  and  generation 
within  the  cell.  The  merits  and  demerits  of  each  method  are  discussed.  For  a  one-dimensional  case  study,  both  methods  were  compared 
quantitatively  and  results  show  that  both  models  accurately  predict  the  polarization  effects  and  water  management  requirements. 

©  2005  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

It  is  believed  that  there  will  be  a  time  in  the  future  when 
global  energy  demands  will  be  met  by  some  source  other  than 
fossil  fuels.  It  is  believed  that  hydrogen  will  play  a  major  role 
in  such  a  future  [1].  The  concept  of  a  hydrogen  economy  de¬ 
scribes  an  economy  where  the  principal  source  of  energy  is 
hydrogen  related.  Fuel  cells,  in  particular  proton  exchange 
membrane  fuel  cells  (PEMFC),  are  expected  to  play  a  major 
role  in  a  future  hydrogen  economy.  Fuel  cells  are  particularly 
attractive  for  use  in  vehicles  as  a  replacement  to  the  com¬ 
bustion  engine.  The  low  temperature  operation  of  a  PEMFC 
(typically  <90  °C)  allows  for  easy  start  up  and  quick  response 
to  changes  in  load  and  operating  conditions. 

However,  a  number  of  issues  need  to  be  resolved  before 
fuel  cells  can  be  commercially  viable.  Typical  proton  ex¬ 
change  membranes  require  precise  water  management,  which 
is  difficult  under  the  variable  load  associated  with  automobile 
driving  conditions.  Dehydration  of  the  membrane  results  in 
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lower  ionic  conductivity  as  well  as  the  risk  of  de-adhesion 
of  the  membrane,  whereas  excessive  water  production  (at 
high  current  densities)  results  in  mass  transport  limitations 
on  the  cathode  side.  Sluggish  electrode  kinetics  also  poses 
a  problem.  The  rate  of  oxygen  reduction  at  the  cathode 
is  much  slower  than  hydrogen  oxidation  at  the  anode,  and 
this  limits  the  performance  of  the  cell.  Also,  trace  amounts 
of  carbon  monoxide  in  the  hydrogen  feed  have  a  deleteri¬ 
ous  effect  on  the  platinum  based  catalyst  typically  used  in 
PEMFCs. 

Fuel  cell  modeling  has  received  much  attention  over  the 
past  15  years  in  an  attempt  to  better  understand  the  phenom¬ 
ena  occurring  within  the  cell.  Parametric  models  allow  engi¬ 
neers  and  designers  to  predict  the  performance  of  the  fuel  cell 
given  geometric  parameters,  material  properties  and  operat¬ 
ing  conditions,  such  as  temperature,  pressure  and  humidity. 
Such  models  are  advantageous  because  experimentation  is 
costly  and  time  consuming.  Furthermore,  experimentation  is 
limited  to  designs,  which  already  exist,  thus  does  not  facili¬ 
tate  innovative  design.  Given  the  highly  reactive  environment 
within  the  fuel  cell,  it  is  often  impossible  to  measure  critical 
parameters,  such  as  temperature,  pressure  and  potential  gra¬ 
dients,  or  species  concentration  within  the  cell.  Thus,  detailed 
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Nomenclature 

a 

effective  catalyst  surface  area  per  unit  volume 
(m-1) 

c 

concentration  (molm-3) 

Cl,  Cl 

integration  constants 

Di 

diffusion  constant  of  species  i  (m2  s-1) 

DiJ 

diffusivity  of  gas  pair  i—j  in  a  mixture  (m2  s  - 1 ) 

E 

potential  (V) 

f 

F/RT(V~l ) 

F 

Faraday’s  constant  (C  equivalent-1) 

i 

ionic  current  density  (Am-2) 

I 

cell  current  density  (Am-2) 

kp 

hydraulic  permeability  (m2) 

k(p 

electro-kinetic  permeability  (m2) 

K, 

Henry’s  law  constant  for  species  i 
(Pam3  mol-1) 

P 

pressure  (Pa) 

R 

universal  gas  constant  (J  mol-1  K-1) 

T 

temperature  (K) 

V 

pore  water  velocity  (ms-1) 

Vs 

superficial  water  velocity  (ms-1) 

Xi 

mole  fraction  of  species  i 

Zi 

charge  number  of  species  i 

Greek  letters 

0£ a?  0ic 

anodic  and  cathodic  transfer  coefficients 

y 

concentration  parameter  in  Butler- Volmer 
equation 

8i,k 

volume  fraction  of  species  i  in  region  k 

? 

stoichiometric  flow  ratio 

K 

ionic  conductivity  of  the  membrane 
(mho  m1) 

P 

pore  water  dynamic  viscosity  (kgm-1  s-1) 

P 

water  molar  density  (molm-3) 

O 

electrical  conductivity  (mhom-1) 

<P 

potential  (V) 

transport  models,  which  accurately  predict  the  flux  and  con¬ 
centration  of  multiple  species  are  required.  Such  information 
is  useful,  for  example,  in  the  loading  of  catalysts.  Transport 
models  can  be  used  to  predict  the  pH  within  the  cell  in  or¬ 
der  to  identify  the  optimum  operating  conditions  for  certain 
catalysts,  and  also  to  identify  where  most  of  the  electrochem¬ 
ical  reactions  take  place.  The  performance  of  ceramic-based 
catalysts  is  pH  dependent.  Transport  models  can  determine 
the  pH  in  the  catalyst  regions  based  on  the  H+  concentration. 
This  would  help  designers  to  optimize  the  cell  for  effective 
catalyst  usage  and  utilization. 

This  paper  reviews  some  of  the  work  done  in  PEMFC 
modeling  over  the  past  15  years,  discusses  contemporary 
trends  and  compares  various  approaches  to  modeling  in  re¬ 
cent  times. 


2.  Categories  of  fuel  cell  models 

A  fuel  cell  model  may  fall  into  one  of  three  cate¬ 
gories:  analytic,  semi-empirical,  or  mechanistic  (theoreti¬ 
cal).  Table  1  categorizes  the  models  reviewed  in  this  paper 
according  to  their  areas  of  investigation  and  dimension  of 
study. 

2.7.  Analytical  models 

Examples  of  analytical  modeling  are  those  reported 
by  Standaert  et  al.  [2,3].  Many  simplifying  assumptions 
were  made  concerning  variable  profiles  within  the  cell  in 
order  to  develop  an  approximate  analytical  voltage  versus 
current  density  relationship.  This  model  also  predicted  water 
management  requirements.  This  was  done  in  the  case  of 
isothermal  and  non-isothermal  cells.  However,  analytical 
models  are  only  approximate  and  do  not  give  an  accurate 
picture  of  transport  processes  occurring  within  the  cell. 
They  are  limited  to  predicting  voltage  losses  and  water 
management  requirements  for  simple  designs.  They  may 
be  useful  if  quick  calculations  are  required  for  simple 
systems. 

2.2.  Semi-empirical  models 

Semi-empirical  modeling  combines  theoretically  derived 
differential  and  algebraic  equations  with  empirically  deter¬ 
mined  relationships.  Empirical  relationships  are  employed 
when  the  physical  phenomena  are  difficult  to  model  or  the 
theory  governing  the  phenomena  is  not  well  understood. 
Springer  et  al.  [4]  developed  a  semi-empirical  model  for 
use  in  a  fuel  cell  with  a  partially  hydrated  membrane  (as 
opposed  to  a  fully  hydrated  membrane).  Empirically  deter¬ 
mined  relationships  were  developed  correlating  membrane 
conductivity  and  electrode  porosity  with  water  content  in  the 
Nation®  membrane.  Most  of  the  models  subsequently  devel¬ 
oped  used  these  correlations  to  determine  the  conductivity  of 
the  Nation®  membrane. 

Amphlett  et  al.  [5]  used  semi-empirical  relationships  to 
estimate  the  potential  losses  and  to  fit  coefficients  in  a  for¬ 
mula  used  to  predict  the  cell  voltage  given  the  operating  cur¬ 
rent  density.  This  model  accounted  for  activation  and  ohmic 
overpotentials.  The  partial  pressures  and  dissolved  concen¬ 
trations  of  hydrogen  and  oxygen  were  determined  empirically 
as  a  function  of  temperature,  current  density  and  gas  chan¬ 
nel  mole  fractions.  Subsequently,  the  reversible  cell  voltage, 
activation  overpotentials  and  cell  resistance  were  correlated 
with  temperature,  partial  pressures,  dissolved  concentrations 
and  operating  current  density.  Pisani  et  al.  [6]  also  used  a 
semi-empirical  approach  to  study  the  activation  and  ohmic 
losses  as  well  as  transport  limitations  at  the  cathode  reactive 
region. 

Maggio  et  al.  [7]  studied  the  water  transport  in  a  fuel 
cell  using  a  semi-empirical  approach.  They  modeled  the 
concentration  overpotential  effect  by  allowing  the  cathode 
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Table  1 


PEMFC  model  categorization  based  on  areas  of  investigation 


Feature 

Analytical 

Semi-empirical 

Mechanistic 

2 

3 

4 

5 

6 

1 

8 

9 

10 

11 

12 

14 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

Dimension 

1 

1 

2 

3 

3 

1 

2 

3 

3 

3 

3 

1 

2 

2 

3 

3 

3 

2 

3 

3 

3 

Polarization 

■j 

V 

■j 

sj 

sj 

■J 

V 

V 

v 

V 

V 

V 

V 

•j 

v 

V 

V 

V 

V 

V 

•j 

v 

V 

V 

V 

v 

V 

V 

V 

Transport  phenomena 

V 

V 

V 

V 

V 

V 

V 

V 

V 

V 

V 

V 

V 

V 

V 

V 

V 

V 

V 

Thermal  effects 

V 

V 

V 

V 

Water  management 

V 

V 

V 

v 

Concentration  effects 

V 

V 

V 

CO  kinetics 

V 

V 

Catalyst  utilization 

V 

V 

Flow  field  effects 

sj 

V 

V 

V 

sj 

V 

V 

■j 

Membrane  conductivity 

•j 

FC  stacks 

V 

V 

gas  porosity  to  be  an  empirical  function  of  current  den¬ 
sity  (since  current  density  is  related  to  water  production). 
The  effective  gas  porosity  was  assumed  to  decrease  lin¬ 
early  with  increasing  current  density.  This  is  due  to  the  in¬ 
creasing  percentage  of  gas  pores  occupied  by  liquid  water. 
Their  results  indicate  that  dehydration  of  the  membrane  is 
likely  to  occur  on  the  anode  side  rather  than  the  cathode 
side. 

Chan  et  al.  [8]  studied  the  effect  of  CO  kinetics  in  the 
hydrogen  feed  on  the  anode  reactive  region.  When  hydro¬ 
gen  is  obtained  from  reformed  fuel,  there  are  trace  amounts 
of  CO  present  which  act  as  poison  to  the  platinum  catalyst. 
The  CO  is  preferentially  adsorbed  onto  the  catalyst  surface 
instead  of  hydrogen,  whereby  decreasing  the  catalyst  surface 
area  available  for  hydrogen  dissociation.  An  empirical  fac¬ 
tor  was  determined  which  represented  the  fraction  of  catalyst 
sites  occupied  by  CO  at  the  anode.  The  result  is  larger  acti¬ 
vation  overpotentials  on  the  anode  side  due  to  slow  electrode 
kinetics. 

Semi-empirical  modeling  has  also  been  used  to  model 
fuel  cell  stacks.  Maxoulis  et  al.  [9]  used  such  an  approach 
to  model  a  fuel  cell  stack  during  automobile  driving  cy¬ 
cles.  They  combined  the  model  of  Amphlett  et  al.  [5]  with 
the  commercial  software  ADVISOR,  which  was  used  to 
simulate  vehicular  driving  conditions.  They  studied  the  ef¬ 
fects  of  the  number  of  cells  per  stack,  electrode  kinetics 
and  water  concentration  in  the  membrane  on  the  fuel  con¬ 
sumption.  They  concluded  that  a  larger  number  of  cells  per 
stack  result  in  greater  stack  efficiency  resulting  in  better  fuel 
economy. 

Semi-empirical  models  are,  however,  limited  to  a  narrow 
corridor  of  operating  conditions.  They  cannot  accurately  pre¬ 
dict  performance  outside  of  that  range.  They  are  very  useful 
for  making  quick  predictions  for  designs  that  already  exists. 
They  cannot  be  used  to  predict  the  performance  of  innovative 
designs,  or  the  response  of  the  fuel  cell  to  parameter  changes 
outside  of  the  conditions  under  which  the  empirical  relation¬ 
ships  were  developed.  Empirical  relationships  also  do  not 
provide  an  adequate  physical  understanding  of  the  phenom¬ 
ena  inside  the  cell.  They  only  correlate  output  with  input. 


2.3.  Mechanistic  models 

Mechanistic  modeling  has  received  the  most  atten¬ 
tion  in  the  literature.  In  mechanistic  modeling,  differ¬ 
ential  and  algebraic  equations  are  derived  based  on  the 
physics  and  electro-chemistry  governing  the  phenomena 
internal  to  the  cell.  These  equations  are  solved  using 
some  sort  of  computational  method.  Mechanistic  models 
can  be  subcategorized  as  multi-domain  models  or  single¬ 
domain  (or  unified)  models.  Fig.  1  gives  a  chronology 
of  the  development  of  mechanistic  modeling.  It  shows 
the  evolution  of  PEMFC  modeling  as  it  increased  in 
complexity. 

2.3.1.  Multi- domain  approach 

Multi-domain  models  involve  the  derivation  of  different 
sets  of  equations  for  each  region  of  the  fuel  cell,  namely  the 
anode  and  cathode  gas  diffusion  regions,  anode  and  cath¬ 
ode  gas  flow  channels,  membrane  and  catalyst  layers.  These 
equations  are  solved  separately  and  simultaneously. 

One  of  the  early  mechanistic  models  for  a  PEMFC  was  the 
pioneering  work  of  Bernardi  and  Verbrugge  [10,11].  They 
developed  a  one-dimensional,  steady  state,  isothermal  model 
which  described  water  transport,  reactant  species  transport, 
as  well  as  ohmic  and  activation  overpotentials.  Their  model 
assumed  a  fully  hydrated  membrane  at  all  times,  and  thus  cal¬ 
culated  the  water  input  and  removal  requirements  to  maintain 
full  hydration  of  the  membrane.  The  model  equations  were 
derived  using  the  Stefan  Maxwell  equations  to  describe  gas 
phase  diffusion  in  the  electrode  regions,  the  Nernst-Planck 
equation  to  describe  dissolved  species  fluxes  in  the  membrane 
and  catalyst  layers,  the  Butler  Volmer  equation  to  describe 
electrode  rate  kinetics  and  Schlogl’s  equation  for  liquid  wa¬ 
ter  transport. 

This  model  was  used  primarily  to  predict  the  polarization 
effects  (due  to  ohmic  and  activation  overpotentials)  and  the 
water  management  requirements.  The  model  computed  the 
required  water  input  at  the  anode  side  and  required  water 
removal  rate  at  the  cathode  side  necessary  to  maintain  full 
hydration  of  the  Nation®  membrane  at  all  times. 
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27 
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19 
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21 
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25 
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18 
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13 

14 
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17 

2D 

24 

12 
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23 

ID 
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10 

11 
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1990  1995  2000  2005 
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Fig.  1.  PEMFC  mechanistic  modeling  evolution. 


Their  model  also  predicted  the  dissolved  hydrogen  and 
oxygen  concentrations  within  the  catalyst  layers  and  it  was 
found  that  at  sufficiently  high  current  densities  most  of  the 
electrochemical  reactions  occurred  on  the  outer  surface  of  the 
catalyst  layer.  This  information  is  vital  to  designers  of  fuel 
cells.  It  allows  them  to  economically  distribute  the  catalyst 
where  it  is  most  needed.  Considering  that  the  platinum  cata¬ 
lyst  is  one  of  the  largest  expenses  in  a  fuel  cell,  this  could  help 
reduce  the  cost.  Although  their  model  was  basic  and  many 
improvements  have  been  made  since,  this  work  served  as  a 
foundation  for  PEMFC  modeling. 

Gurau  et  al.  [12]  developed  a  two-dimensional  model  to 
determine  species  concentrations  within  the  fuel  cell  and  the 
effect  on  fuel  cell  performance  of  gas  diffuser  porosity,  air 
flow  rates  in  the  gas  channel  and  temperature.  This  model 
is  based  on  the  multi-domain  approach  with  three  domains 
considered:  gas  diffusers/gas  flow  channels,  catalyst  layers 
and  membranes.  The  gas  diffusers  and  gas  flow  channels 
were  combined  into  one  domain  by  writing  the  governing 
equations  for  each  region  in  a  similar  form.  As  a  result  the 
solution  methodology  was  able  to  accommodate  both  regions 
into  one  domain. 

2.3.2.  Single -domain  approach 

Gurau  et  al.  [12]  showed  that  since  the  governing  differen¬ 
tial  equations  in  the  gas  flow  channels  and  the  gas  diffusion 
electrodes  are  similar,  the  equations  can  be  combined  for  both 
regions.  The  computational  effect  is  that  both  regions  can  be 
considered  as  one  domain  where  no  internal  boundary  con¬ 


ditions  or  statements  of  continuity  need  be  defined.  The  only 
difference  is  that  material  properties  and  source  terms  assume 
different  values  for  the  two  regions.  This  forms  the  basis  of 
the  single-domain  approach. 

Instead  of  combining  two  regions  into  one  domain,  the 
single-domain  approach  combines  all  the  regions  of  interest 
into  one  domain.  Conservation  equations  are  defined  which 
govern  the  entire  domain  of  interest,  typically  the  entire  fuel 
cell  (gas  flow  regions  and  the  membrane  electrode  assembly). 
In  each  region,  the  differences  are  accounted  for  by  source  and 
sink  terms.  All  equations  are  written  in  the  form  of  a  generic 
convection-diffusion  equation,  and  all  terms,  which  do  not 
fit  that  format  are  dumped  into  the  source  or  sink  term.  This 
formulation  allows  for  solution  using  known  computational 
fluid  dynamics  (CFD)  methods. 

The  principles  of  CFD  were  first  applied  to  fuel  cells 
by  Wang  et  al.  [13]  and  Zhou  and  Fiu  [14].  In  the  unified 
approach,  all  the  governing  differential  equations  were  ar¬ 
ranged  into  a  standard  form,  which  could  be  discretized  using 
the  principles  of  CFD  [15]  or  solved  using  a  CFD  software 
package. 

d& 

—  =  V(n0)  +  V(CV0)  +  5 
dt 

The  respective  terms  in  above  equation  represent  transient, 
convection,  diffusion  and  source.  The  general  variable  may 
refer  to  potential,  temperature,  pressure,  velocity,  concentra¬ 
tion  or  phase  fraction.  For  steady  state  operation  the  transient 
term  vanishes.  This  general  equation  is  a  conservation  equa¬ 
tion,  where  the  source  term  represents  material  consumption 
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and  generation  (at  catalyst  layers),  phase  change  and  any 
other  term,  which  cannot  fit  the  general  format,  but  must  be 
accounted  for.  The  application  of  CFD  to  fuel  cells  has  paved 
the  way  for  the  subsequent  development  of  multidimensional 
models. 


3.  Modeling  considerations 

Mechanistic  modeling  (single  and  multi-domain)  has 
been  utilized  to  study  a  wide  range  of  phenomena 
including  polarization  effects  (activation,  ohmic  and  con¬ 
centration  overpotentials),  water  management,  thermal  man¬ 
agement,  CO  kinetics,  catalyst  utilization  and  flow  field 
geometry. 

3.1.  Parametric  models 

All  models  are  parametric  in  that  they  predict  the  output 
performance  for  various  input  parameters,  typically  temper¬ 
ature,  pressure  and  humidity.  Wang  et  al.  [16]  developed  a 
three-dimensional  parametric  model,  considering  the  effects 
of  temperature,  humidity  and  pressure.  It  was  found  that  the 
performance  of  the  fuel  cell  improved  with  increasing  tem¬ 
perature  if  the  inlet  gases  are  fully  humidified.  If  the  gases  are 
not  fully  humidified,  dehydration  of  the  membrane  is  likely 
to  occur  resulting  in  reduced  conductivity  values,  hence  re¬ 
duced  cell  performance.  They  also  found  that  at  low  current 
densities  anode  humidification  is  required,  but  not  at  higher 
current  densities.  This  is  because  at  high  current  densities, 
sufficient  water  is  produced  at  the  cathode  to  keep  the  mem¬ 
brane  hydrated. 

Their  results  further  show  that  cathode  humidification  is 
not  significant  at  all,  especially  at  high  current  densities.  This 
is  because  dehydration  is  likely  to  occur  on  the  anode  side 
and  flooding  on  the  cathode  side.  Therefore,  humidifying  the 
cathode  gas  stream  adds  no  benefit. 

Finally,  increasing  the  pressure  of  the  inlet  gases  was  seen 
to  improve  performance  by  increasing  the  activation  currents 
and  the  partial  pressures  of  reactant  gases. 

The  authors  report  that  at  higher  current  densities,  their 
model  overestimates  the  cell  current  density  compared  to  ex¬ 
perimental  results.  The  reason  for  this  is  that  the  model  did 
not  take  into  account  mass  transport  effects. 

3.2.  Mass  transport  effects 

Mass  transport  limitations  or  concentration  overpoten¬ 
tials  are  caused  when  the  reactants  cannot  be  supplied  fast 
enough  for  the  required  rate  of  chemical  reaction  to  take 
place.  This  happens  especially  at  high  current  densities  when 
large  amounts  of  liquid  water  are  produced  at  the  cathode. 
Liquid  water  has  a  two-fold  effect.  It  dilutes  the  reactants 
thus  reducing  its  concentration  near  the  catalyst  sites,  and  it 
reduces  the  effective  gas  porosity  thus  “blocking”  the  reac¬ 
tants  from  reaching  the  catalyst  layer. 


Baschuk  and  Li  [17]  developed  a  one-dimensional  model, 
which  accounted  for  cathode  mass  limitation  effects  by  al¬ 
lowing  variable  degrees  of  flooding  at  the  cathode  catalyst 
layer/backing  region.  They  account  for  concentration  over¬ 
potential  as  a  result  of  the  decreased  concentration  of  dis¬ 
solved  oxygen  in  the  catalyst  region  due  to  the  excessive 
water  content.  Darcy’s  law  is  used  to  obtain  the  drop  in  par¬ 
tial  pressure  of  the  oxygen  at  the  cathode  catalyst  layer,  and 
Henry’s  law  is  used  to  determine  the  dissolved  oxygen  con¬ 
centration.  Their  modeled  results  showed  excellent  agree¬ 
ment  with  experimental  results.  The  model  also  predicted  that 
increasing  the  cell  pressure  lowers  the  limiting  current  den¬ 
sity.  High  pressures  result  in  maximum  flooding  occurring 
at  lower  current  densities,  and  this  effect  is  more  significant 
than  the  increase  in  partial  pressure  of  the  oxygen.  The  results 
also  showed,  predictably,  that  increasing  the  temperature  in¬ 
creases  the  limiting  current  density. 

Um  et  al.  [18]  developed  a  transient  model  based  on  the 
unified  approach,  which  studied  the  effects  of  hydrogen  dilu¬ 
tion  along  the  anode  gas  channel.  The  two-dimensional  model 
considered  flow  perpendicular  to  the  membrane  electrode  as¬ 
sembly  (MEA)  cross-section,  as  well  as  in  the  direction  of 
flow  in  the  gas  channels.  As  hydrogen  diffuses  from  the  gas 
channel  into  the  gas  diffusion  region,  its  concentration  along 
the  gas  channel  decreases  resulting  in  a  two-dimensional  con¬ 
centration  gradient  in  the  gas  diffusion  electrodes.  The  result 
is  that  mass  transport  limitations  are  seen  on  the  anode  side 
especially  at  high  current  densities,  and  when  reformed  fuel 
is  used  instead  of  pure  hydrogen.  At  high  current  densities, 
hydrogen  is  extracted  from  the  flow  channels  at  a  much  fast 
rate  than  at  low  current  densities.  With  reformed  fuel,  the 
partial  pressure  of  the  hydrogen  is  already  lowered  by  the 
presence  of  carbon  dioxide  in  the  gas  feed,  so  as  it  is  used  up 
toward  the  end  of  the  gas  channel,  the  partial  pressure  of  the 
hydrogen  may  be  too  low  and  it  may  not  be  able  to  diffuse 
to  the  anode  catalyst  layer  fast  enough.  The  result  is  anode 
side  mass  transport  limitations.  Such  phenomena  cannot  be 
studied  using  one-dimensional  models. 

Zhou  and  Liu  [19]  developed  a  3-D  model  of  a  PEMFC 
taking  into  account  CO  effects  in  the  anode  gas  stream.  The 
model  accounts  for  poisoning  of  the  catalyst  as  well  as  hy¬ 
drogen  dilution  due  to  the  inert  gases.  An  interesting  result  is 
that  along  the  anode  gas  channel  the  hydrogen  concentration 
increases  in  the  direction  of  flow  because  the  CO  is  depleted 
at  a  faster  rate  due  to  preferential  adsorption  at  the  catalyst 
sites.  They  also  found  that  the  optimum  porosity  of  the  gas 
diffusion  layer  is  much  lower  for  a  fuel  cell  using  reformate 
than  one  using  pure  hydrogen. 

3.3.  Thermal  management 

The  electrochemical  reactions  taking  place  in  a  fuel  cell 
are  exothermic,  i.e.,  they  give  off  heat.  Heat  is  also  produced 
by  irreversibilities  in  the  cell  such  as  activation  losses  and 
ohmic  effects.  Heat  removal  is  a  critical  design  issue  for  fuel 
cells.  Excessive  heat  generation  may  result  in  dehydration  of 
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the  membrane  resulting  in  decreased  conductivity,  and  ther¬ 
mal  stresses  resulting  in  mechanical  failure.  Various  models 
have  accounted  for  non-isothermal  effects  and  heat  transfers. 

Berning  et  al.  [20,21]  used  the  unified  approach  to  de¬ 
velop  a  three-dimensional  non-isothermal  fuel  cell  model. 
The  model  studied  reactant  concentrations,  current  density 
distributions  and  temperature  gradients  within  the  cell  as  well 
as  water  flux  and  species  transport.  For  gas  flow  fields  sepa¬ 
rated  by  current  collecting  plates,  three-dimensional  effects 
were  observed  due  to  the  unevenness  of  the  hydrogen  and 
oxygen  supply.  These  effects  were  pronounced  under  the  col¬ 
lector  plate  land  areas.  These  effects  may  result  in  transport 
limiting  conditions  at  high  current  densities. 

The  development  of  a  non-isothermal  model  was  intended 
to  study  the  heat  transfers  within  the  cell.  It  was  observed  that 
a  temperature  difference  of  2-3  K  existed  within  the  cell.  Yan 
et  al.  [22]  performed  a  similar  study  and  found  a  temperature 
variation  within  the  cell  of  the  same  magnitude.  However,  the 
magnitude  of  the  heat  transfer  was  not  reported  so  it  is  difficult 
to  compare  the  magnitude  of  the  conductive  heat  transfer 
relative  to  the  total  heat  transfer.  This  information  would  have 
helped  justify  the  need  for  non-isothermal  modeling.  If  the 
heat  transfer  by  conduction  were  small  compared  to  other 
heat  transfers,  then  a  temperature  difference  of  2-3  K  could 
hardly  be  significant. 

Wohr  et  al.  [23]  investigated  heat  management  for  fuel  cell 
stacks.  In  their  one-dimensional,  non-isothermal  model  they 
considered  the  gas  diffusion  region  to  be  a  homogeneous  dis¬ 
tribution  of  cylindrical  pores,  through  which  transport  was 
governed  by  the  “dusty  gas  model”.  Water  transport  in  the 
electrodes  was  assumed  to  occur  by  surface  diffusion  or  cap¬ 
illary  effect.  They  considered  the  heating  effects  due  to  ohmic 
resistance  in  the  membrane  and  heat  generation  due  to  the  en¬ 
tropy  of  reaction.  In  the  absence  of  any  heat  removal  strategy, 
the  temperature  difference  in  a  stack  of  four  cells  was  8  K. 
The  temperatures  of  the  innermost  cells  in  the  stack  were  the 
highest.  The  effect  is  that  the  membranes  of  the  inner  cells 
would  dehydrate. 

3.4.  Water  management 

Water  management  is  another  critical  aspect  of  PEM  fuel 
cells.  At  high  current  densities,  excessive  water  transport 
across  the  membrane  and  water  production  at  the  cathode 
result  in  flooding  of  the  electrodes  and  mass  transport  limita¬ 
tions.  At  low  current  densities,  dehydration  of  the  membrane 
may  occur  at  the  anode  side.  Water  must  be  supplied  to  the 
fuel  cell  at  the  anode  and  removed  at  the  cathode  in  order  to 
maintain  effective  membrane  humidification. 

The  Bernardi  and  Verbrugge  model  [10,11]  determines 
the  required  water  addition  and  removal  necessary  to 
maintain  full  membrane  humidification  at  a  given  current 
density.  Fuller  and  Newman  [24]  developed  a  pseudo  two- 
dimensional  model,  which  predicted  water  and  thermal  man¬ 
agement  as  well  as  fuel  utilization  for  a  fuel  cell  operating 
with  reformed  methanol  as  the  fuel. 


The  model  of  Wohr  et  al.  [23]  showed  that  for  fuel  cell 
stacks  water  management  becomes  even  more  difficult  and 
is  strongly  related  to  thermal  management.  The  temperatures 
of  the  inner  cells  of  the  stack  are  higher  than  the  outer  cells 
resulting  in  membrane  dehydration.  Water  management  is 
strongly  interrelated  with  thermal  management/heat  removal. 
Other  strategies  for  effective  water  management  involve  the 
geometry  of  the  gas  flow  field,  e.g.,  counter  flow  versus  co¬ 
flow,  and  using  serpentine  and  interdigitated  flow  fields  rather 
than  straight  channel  flow  fields. 

3.5.  Flow  field  geometry 

One  of  the  advantages  of  computational  modeling  over 
experimentation  to  assess  the  performance  of  a  fuel  cell  is  the 
ability  to  evaluate  innovative  designs.  One  area  where  this  is 
evident  is  in  the  consideration  of  flow  field  geometry.  Models 
have  been  developed  for  straight  flow  channels,  serpentine 
flow  channels  and  interdigitated  flow  fields. 

Ge  and  Yi  [25]  developed  a  two-dimensional  model  to 
study  the  effects  of  flow  mode  in  straight  gas  channels,  i.e., 
counter  flow  versus  co-flow.  It  was  found  that  the  flow  mode 
only  made  a  difference  when  dry  or  low  humidity  inlet  gases 
were  used.  For  such  cases,  counter  flow  operation  produced 
better  results  since  by  so  doing  the  reactant  gases  were  suf¬ 
ficiently  humidified  internally.  If  the  inlet  gases  are  already 
humidified,  the  flow  mode  makes  little  difference.  The  reason 
for  this  is  that  for  high  humidity  gases,  the  increase  in  mem¬ 
brane  conductivity  due  to  the  high  humidity  is  counteracted 
by  the  increase  in  cathode  concentration  overpotential  due  to 
the  presence  of  liquid  water.  This  is  the  case  whatever  the  flow 
mode.  However,  for  low  humidity  gases,  counter  flow  opera¬ 
tion  allows  for  internal  humidification  of  the  gas  streams.  For 
co-flow  low  humidity  gases,  the  membrane  dehydrates.  This 
information  gives  the  designers  of  fuel  cells  an  alternative  to 
humidifying  the  gas  streams. 

Dutta  et  al.  [26]  used  the  unified  approach  to  study  mass 
transport  between  the  channels  of  a  PEMFC  with  a  serpen¬ 
tine  flow  field.  Their  model  is  three-dimensional  and  allows 
for  multi-species  transport.  They  studied  the  effect  of  flow 
channel  width  in  the  serpentine  flow  field  on  velocity  distri¬ 
bution,  gas  mixture  distribution  and  reactant  consumption. 
Serpentine  flow  fields  allow  for  a  greater  area  for  diffusion 
of  the  supply  gases.  Their  results  show  that  for  low  humidity 
conditions,  water  transport  is  dominated  by  electro- osmotic 
effects,  i.e.,  water  flows  from  anode  to  cathode  at  the  side  of 
the  cell  closer  to  the  gas  channel  inlet.  At  the  outlet  side  of 
the  cell,  water  transport  is  dominated  by  back  diffusion,  and 
it  flows  in  the  opposite  direction.  Thus  the  serpentine  flow 
field  allows  for  circulation  of  the  water  within  the  cell. 

Nguyen  et  al.  [27]  developed  a  three-dimensional  model 
which  accounts  for  mass  and  heat  transfer,  current  and  po¬ 
tential  distribution  within  a  cell  using  a  serpentine  flow  field. 
Their  results  show  that  oxygen  concentration  along  the  gas 
channels  decrease  in  the  direction  of  flow.  Also,  in  the  gas 
diffusion  layer,  the  oxygen  concentration  is  a  minimum  under 
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the  land  area.  At  high  current  densities  the  oxygen  is  almost 
completely  depleted  under  the  land  areas.  The  result  is  an  un¬ 
even  distribution  of  oxygen  concentration  along  the  catalyst 
layer  resulting  in  local  overpotentials,  which  vary  spatially. 
A  unique  feature  of  this  model  is  a  voltage-to-current  (VTC) 
algorithm,  which  allows  for  the  solution  of  the  potential  field 
and  the  local  activation  overpotential.  Since  the  reactant  con¬ 
centration  is  not  constant  across  the  catalyst  layer,  the  acti¬ 
vation  overpotential  will  not  be  constant.  Their  simulations 
show  a  variation  in  local  activation  overpotential  from  0.31 
to  0.37  V  at  a  current  density  of  1.2  A  cm-2.  This  VTC  al¬ 
gorithm  however,  comes  with  a  computational  cost.  It  slows 
down  the  solution  requiring  6000-8000  iterations  for  conver¬ 
gence. 

Um  and  Wang  [28]  used  a  three-dimensional  model  to 
study  the  effects  an  interdigitated  flow  field.  The  model  ac¬ 
counted  for  mass  transport,  electrochemical  kinetics,  species 
profiles  and  current  density  distribution  within  the  cell.  In¬ 
terdigitated  flow  fields  result  in  forced  convection  of  gases, 
which  aids  in  liquid  water  removal  at  the  cathode.  This 
would  help  improve  performance  at  high  current  densities 
when  transport  limitations  due  to  excessive  water  produc¬ 
tion  are  expected.  The  model  shows  that  there  is  little  to  no 
difference  at  low  to  medium  current  densities  between  an 
interdigitated  flow  field  and  a  conventional  flow  field.  How¬ 
ever,  at  higher  current  densities,  a  fuel  cell  with  an  inter¬ 
digitated  flow  field  has  a  limiting  current,  which  is  nearly 
50%  greater  than  an  equivalent  cell  with  a  conventional  flow 
field.  Because  of  the  flow  field,  three-dimensional  effects  un¬ 
der  the  current  collector  land  area,  known  as  rib  effects,  are 
prominent. 

Seigel  et  al.  [29]  also  modeled  a  fuel  cell  with  an  interdig¬ 
itated  flow  field.  Theirs  was  a  two-dimensional  steady  state 
model,  which  studied  transport  limitations  due  to  water  build 
up  in  the  cathode  catalyst  region.  They  considered  water  in 
three  phases:  liquid,  gas  and  dissolved  (membrane  phase). 
They  found  that  treating  the  catalyst  layer  as  a  very  thin  in¬ 
terface  underestimates  the  transport  limitations  due  to  water 
build-up.  Hence,  they  modeled  the  catalyst  layer  as  a  finite 
region.  Their  model  showed  that  20-40%  of  the  water  build¬ 
ing  up  at  the  cathode  catalyst  layer  comes  from  water  which 
is  transported  across  the  membrane.  This  problem  may  be 
counteracted  by  applying  a  pressure  differential  to  force  back 
diffusion  of  water,  i.e.,  from  cathode  to  anode. 

Using  the  multi-domain  approach  Hu  et  al.  [30,31]  devel¬ 
oped  a  three-dimensional  two-phase  model  for  a  fuel  cell. 
They  gave  boundary  conditions,  which  could  be  used  for 
straight  flow  channels  as  well  as  interdigitated  flow  fields. 
Unlike  previous  models,  which  assume  separate  flow  chan¬ 
nels  for  gases  and  liquids,  this  model  assumes  a  two-phase 
mixture.  Water  properties  such  as  specific  volume  change 
depending  on  the  degree  of  mixture.  They  used  a  CFD  algo¬ 
rithm  to  solve  for  the  flow  field  in  the  gas  flow  channels  and 
diffusion  regions,  and  the  fourth  order  Runge-Kutta  method 
together  with  a  shooting  technique  to  solve  for  the  flow  field 
in  the  catalyst  layers  and  the  membrane. 


For  the  interdigitated  flow  field,  results  show  that  the  oxy¬ 
gen  concentration  is  higher  and  liquid  water  saturation  is 
lower  than  those  for  a  conventional  straight  channel  flow 
field.  The  higher  oxygen  concentration  results  in  fast  reac¬ 
tion  rates  and  the  lower  liquid  water  saturation  results  in 
less  concentration  overpotential.  It  is  also  shown  that  the 
local  current  densities  are  much  more  uniform  with  an  in¬ 
terdigitated  flow  field  than  with  a  conventional  flow  field. 
However,  the  performance  of  a  fuel  cell  with  an  interdigi¬ 
tated  flow  field  is  only  shown  to  be  better  than  that  for  one 
with  a  conventional  flow  field  if  the  inlet  gases  are  well 
humidified.  This  is  because  the  interdigitated  field  aids  in 
water  removal,  but  does  not  aid  in  hydration  of  an  already 
de-hydrated  membrane.  So,  the  internal  gases  need  to  be 
humidified. 

Using  the  unified  approach,  Kumar  and  Reddy  [32]  studied 
the  effects  of  having  metal  foam  in  the  flow  field  of  the  bipolar 
plates.  Their  three-dimensional  steady  state  model  shows  that 
decreasing  the  permeability  of  the  gas  flow  field  improves 
performance.  This  is  because  at  low  flow  field  permeability 
reactant  gases  are  transported  by  forced  convection  rather 
than  diffusion.  Having  many  tiny  gas  channels  results  in  a 
lower  permeability  than  having  few  large  channels.  However, 
due  to  limitations  in  machining  processes,  the  flow  channels 
cannot  be  made  too  small.  Placing  metal  foam  in  the  flow 
field  allows  the  flow  field  permeability  to  be  lowered  without 
resorting  to  precise  machining  processes.  They  found  that 
decreasing  the  permeability  from  10~6tol0_12m2  increases 
the  “average  current  density”  of  their  system  from  5943  to 
8425  A  m“2. 


4.  Comparisons 

In  this  section,  we  compare  the  unified  (single-domain) 
approach  with  the  multi-domain  approach  to  solving  the  gov¬ 
erning  equations.  In  Section  3,  we  saw  that  both  approaches 
have  been  used  to  solve  three-dimensional  problems  with 
very  complex  flow  fields.  The  early  mechanistic  models  were 
all  solved  using  the  multi-domain  approach.  With  the  in¬ 
troduction  of  CFD  methods  to  fuel  cell  modeling,  the  door 
opened  for  multi-dimensional  modeling.  Most  of  the  mod¬ 
els  developed  over  the  past  5  years  were  solved  using  the 
single-domain  approach. 

Although  both  approaches  have  been  used  for  multidi¬ 
mensional  models,  the  single-domain  approach  more  easily 
lends  itself  to  multidimensional  modeling.  With  the  multi- 
domain  approach,  internal  boundary  conditions  or  conditions 
of  continuity  must  be  specified  at  each  interface  between  re¬ 
gions,  which  could  become  cumbersome  in  two-  and  three- 
dimensions. 

The  single-domain  approach  also  more  easily  lends  itself 
to  be  implemented  in  commercial  CFD  codes  since  the  so¬ 
lution  methods  to  CFD  problems  are  well  established.  The 
solution  methods  for  the  multi-domain  equations  are  not  as 
standardized  as  those  for  the  single-domain  equations.  As  a 
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result,  the  model  development  time  is  shortened  using  the 
single-domain  approach. 

An  interesting  observation  is  that  many  of  the  earlier 
models  developed  by  researchers  with  a  chemical  engineer¬ 
ing/chemistry  background  were  based  on  the  multi-domain 
approach,  whereas  most  of  the  newer  models  developed  by 
researchers  with  a  mechanical  engineering  background  are 
based  on  the  unified  approach. 

It  is  difficult  to  say  which  method  converges  faster  since 
that  depends  heavily  on  what  is  being  modeled,  and  how 
efficiently  the  programs  are  written.  It  has  been  reported  that 
the  single-domain  approach  requires  longer  solution  times 
[29],  however,  a  possible  reason  is  that  commercial  codes 
are  very  general  and  not  computationally  optimized  for  any 
specific  situation. 
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5.  Case  study 

For  the  purpose  of  comparison,  we  solved  the  model  devel¬ 
oped  by  Bernardi  and  Verbrugge  [11]  using  both  approaches. 
The  same  terminology  as  employed  by  the  authors  [11]  is 
used  together  with  the  same  material  properties  and  operat¬ 
ing  conditions.  The  model  is  one-dimensional,  the  only  di¬ 
mension  considered  is  that  perpendicular  to  the  MEA  cross 
section.  The  solution  domain  consists  of  the  MEA,  i.e.,  the 
anode  and  cathode  gas  diffusion  regions,  catalyst  layers  and 
the  membrane.  The  assumptions  of  the  model  are: 


I 

— x 

4  F 


1  T  X]Sf2 


/  neff 

/  ^W,Q  2 

l  Deff 

V  ^W,N2 


(6) 

(7) 

(8) 

(9) 


•  constant  cell  temperature  (isothermal),  in  this  case  353  K; 

•  ideal  and  well-mixed  gases  in  the  gas  chambers; 

•  steady  state  operation; 

•  total  gas  pressure  within  each  diffuser  region  is  taken  to 
be  constant  since  gas  phase  viscosity  is  small  compared  to 
liquid  phase  viscosity; 

•  inlet  gas  streams  (hydrogen  and  air)  are  saturated  and  gases 
are  assumed  to  be  fully  saturated  throughout  the  cell; 

•  separate  flow  channels  for  gases  and  liquids  in  the  porous 
regions  of  the  gas  diffusion  layers,  i.e.,  liquid  and  gas 
phases  will  be  treated  as  separate  single  phase  flows  rather 
than  a  two  phase  mixture. 

The  Nernst-Planck  equation  is  used  to  describe  the  trans¬ 
port  of  dissolved  species  in  the  membrane  and  catalyst  layers. 
Transport  of  dissolved  species  occurs  via  diffusion,  convec¬ 
tion  and  migration.  Schlolgl’s  equation  is  used  to  describe  the 
transport  of  liquid  water  via  electro-osmotic  drag  and  back 
diffusion.  The  Stefan  Maxwell  equations  are  used  to  describe 
the  diffusion  of  gases  within  the  gas  diffusion  regions.  The 
Butler- Volmer  equation  is  used  to  describe  electrochemical 
rate  kinetics  in  the  catalyst  layers. 

5.7.  Multi-domain  approach 

The  multi-domain  model  is  characterized  by  equations 
(1)— (25)  in  each  region  as  reported  in  [1 1].  Note  that  in  these 
equations,  the  current  density,  7,  assumes  a  positive  value.  In 
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The  boundary  conditions  or  conditions  of  continuity  must 
be  specified  at  each  interface.  At  the  anode  gas  chamber/gas 
diffuser  interface: 


Vs  =  VO 


P  =  P  0 


(26) 

(27) 


where  vo  represents  the  required  water  input  at  the  anode  in 
order  to  maintain  full  membrane  hydration  at  all  times.  At 
the  anode  diffuser/catalyst  interface: 
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due  to  the  condensation  of  vapor  to  liquid 
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i.e.,  continuous  flux  of  dissolve  hydrogen 
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cq2  =0  i.e. ,  no  dissolved  oxygen  crossover  in  membrane 

(31) 

Similarly  at  the  membrane/cathode  catalyst  layer  interface: 
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At  the  cathode  catalyst  layer/gas  diffuser  interface: 
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The  algorithm  used  to  solve  the  given  system  of  equations 
is  shown  in  Fig.  2.  The  fourth  order  Runge-Kutta  method  is 
used  to  solve  the  equations  in  each  region.  Shooting  tech¬ 
niques  were  used  to  obtain  convergence.  Three  unknown 
terms  had  to  be  initially  guessed:  no,  C\  and  C2.  C\  and  C2  are 
parameters  used  in  the  solution  of  the  second  order  equations 
in  the  membrane  and  catalyst  regions.  So,  three  convergence 
loops  were  required.  The  program  must  converge  to  the  ap¬ 
propriate  value  of  no  such  that  the  pressure  variable  matches 
the  end  boundary  conditions.  Similarly,  the  appropriate  val¬ 
ues  of  Ci  and  C2  must  be  found  for  the  dissolved  hydrogen 
and  oxygen  values  to  match  the  given  internal  boundary  con¬ 
ditions.  A  Newton-Raphson  method  was  used  together  with 
the  shooting  technique  to  speed  up  the  rate  of  convergence. 


GuessVb 

— I  Until  P converges 


Solve  Anode  Gas  Diffuser  Equations 

_ ^ _ 

Solve  Cathode  Gas  Diffuser  Equations 

jk 

Calculate  the  activation  overpotentials 


\|z  GuessCr 


Until  CH2  converges 


Solve  Anode  Catalyst  Layer  Equations 


\|/  Guess  C2 


Until  CO2  converges 


Solve  Cathode  Catalyst  Layer  Equations 


Fig.  2.  Solution  algorithm  for  the  multi-domain  method. 
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5.2.  Single-domain  approach 

The  single-domain  approach  consists  of  conservation 
equations,  which  govern  the  entire  domain.  The  governing 
equations  (38)-(53)  are  shown.  Note  all  units  are  molar  quan¬ 
tities. 

Continuity  equation: 

V(pS)  =  (38) 

where,  Sm  refers  to  the  production  and  consumption  of 
species  at  the  catalyst  layers. 

Sm  2L  =  +  Sw,a  at  the  anode  catalyst  layer  (39) 


>111,0 


=  So ,  +  5. 


w,c 


I 

2F' 


at  the  cathode  catalyst  layer 


i.e. ,  hydrogen  consumed  at  the  anode  catalyst  layer 


(40) 


(41) 


Fig.  3.  Solution  algorithm  for  the  single- domain  method. 

V(pi;.vw)  =  V(7VW)  +  Sw,a  T  $w,c  (31) 

where  the  diffusion  flux  is  given  by  Fick’s  law: 


I 

4 F' 


i.e. ,  oxygen  consumed  at  the  cathode  catalyst  layer 


(42) 


c  l  (  PI  \ 

w’a  2F  \po  ~  Pw  /  ’ 

i.e. ,  water  dissolving  at  the  anode  catalyst  layer  (43) 


I 


i.e.,  water  produced  at  the  cathode  catalyst  layer  (44) 

In  one-dimensional  form,  equations  (41)-(47)  are  identi¬ 
cal  to  equations  (4),  (5),  (8),  (9)  and  (17). 

The  momentum  equation: 


W(pvv  +  P)  +  Sp  —  Sm^v  +  SQ0 


(45) 


where  the  source  term  in  the  momentum  equation  for  porous 
media  is  given  by  Darcy’s  Law: 

SP  =  ~fv  (46) 

/Cp 


The  relatively  small  value  of  hydraulic  permeability  kv  re¬ 
sults  in  the  source  term  dominating  the  momentum  equation 
and  hence  the  problem  reduces  to  equation  (17).  The  electro- 
osmotic  transport  of  liquid  water  in  the  membrane  is  given 
by: 


ZfCfFVcj) 


(47) 


The  transport  of  species  H2,  O2,  N2  and  H2O  vapor  are  given 
by  the  following  conservation  equations  (48)-(51): 


V(/mv  h2)  =  -V(A^h2)  +  5h2  (48) 

V(p2x  o2)  =  —  V(A^o2)  +  So2  (49) 

V(pi5vN2)  =  —  V(A^n2)  (50) 


Ni  =  -pDijWxi  (52) 

Once  again  it  is  obvious  that  these  are  the  same  equations  as 
the  multi-domain  approach  just  written  in  a  different  form. 
Finally  the  potential  field  equation,  which  accounts  for  ohmic 
and  activation  overpotentials,  is  given  by: 

V(a0)  =  -/  +  S<p  (53) 

where  accounts  for  activation  overpotentials  at  both  cata¬ 
lyst  layers. 

These  equations  are  solved  using  a  CFD  finite  volume 
method  as  described  in  Patankar  [15].  The  momentum  and 
continuity  equations  are  solved  first  using  the  SIMPLE  algo¬ 
rithm  with  a  staggered  grid  system.  In  the  pressure  correc¬ 
tion  subroutine,  the  tri-diagonal  matrix  algorithm  (TDMA) 
is  used  to  speed  up  the  solution  process.  Then  the  species 
equations  are  solved,  after  which  the  potential  loss  calcula¬ 
tions  are  made.  The  whole  procedure  is  then  repeated  until 
convergence  is  achieved.  Fig.  3  shows  the  solution  algorithm 
for  the  unified  approach. 


6.  Results  and  discussion 

In  this  section,  results  of  both  models  are  compared.  The 
results  of  interest  are  polarization  curves  and  water  man¬ 
agement  predictions.  Both  models  predict  the  voltage  losses 
hence  the  net  cell  voltage  for  a  given  current  density,  and 
both  predict  the  required  water  input  and  output  required  for 
maintaining  full  hydration  of  the  membrane  at  all  times. 

Fig.  4  shows  the  polarization  curves  obtained  using  both 
models  as  well  as  the  experimental  results  obtained  by  Ti- 
cianelli  et  al.  [33]  for  the  same  operating  conditions.  Saturated 
hydrogen  at  3  atm  was  used  as  fuel  and  saturated  air  (21% 
oxygen,  79%  nitrogen  by  mole  ratio)  at  5  atm  was  used  as 
oxidant.  The  cell  temperature  was  80  °C.  Both  model  curves 
are  virtually  indistinguishable  from  each  other.  This  is  not 
surprising  since  both  models  solve  essentially  the  same  sets 
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Fig.  4.  Polarization  curves,  theoretical  and  experimental. 


of  equations,  just  using  a  different  technique.  Both  compare 
well  with  the  experimental  data  in  the  activation  and  ohmic 
overpotential  regions.  Both  however,  show  signs  of  deviat¬ 
ing  from  the  experimental  results  at  higher  current  densi¬ 
ties.  The  reason  for  this  is  that  the  given  model  does  not 
account  for  concentration  overpotential,  which  is  expected 
to  become  increasingly  prominent  as  the  current  density 
increases. 

Fig.  5  shows  the  water  management  requirements  as  a 
function  of  cell  current  density.  It  shows  both  the  required  in¬ 
put  of  water  at  the  anode  and  the  required  extraction  rate  at  the 
cathode  for  proper  management,  i.e.,  maintaining  full  hydra¬ 


tion  of  the  membrane  at  all  times.  At  lower  current  densities 
liquid  water  flows  in  the  reverse  direction,  i.e.,  out  through 
the  anode  hence  the  negative  values.  The  reason  for  this  as 
explained  in  [1 1]  is  that  at  low  current  densities,  the  effects  of 
back-pressure  are  larger  than  the  effects  of  electro-osmotic 
drag.  Since  the  cathode  gas  pressure  is  higher  than  the  anode 
pressure,  this  pressure  gradient  forces  the  water  to  flow  from 
cathode  to  anode.  At  higher  current  densities,  the  electro- 
osmotic  effects  dominate  and  water  flows  from  the  anode  to 
cathode.  Once  again  both  model  curves  are  very  close.  This 
is  expected  since  the  governing  equations  are  essentially  the 
same. 


Fig.  5.  Water  management  requirements  vs.  current  density. 
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In  this  study,  the  single-domain  solution  converged  much 
faster  than  the  multi-domain  solution.  There  are  a  number  of 
reasons  for  this.  The  finite  volume  method  used  in  the  single¬ 
domain  solution  entailed  a  larger  computational  error  than  the 
single  step  Runge-Kutta  method  used  in  the  multi-domain  so¬ 
lution.  The  Runge-Kutta  method  requires  up  to  four  times  as 
many  calculations.  Secondly,  the  multi-domain  model  treated 
the  catalyst  layer  as  a  finite  sized  region  whereas  the  single¬ 
domain  model  treated  the  catalyst  layers  as  infinitesimally 
small  interfaces  where  species  are  consumed  and  produced. 
For  the  purpose  of  this  study,  the  only  output  parameters 
of  interest  were  the  polarization  curves  and  the  water  man¬ 
agement  information.  Therefore,  the  treatment  of  the  catalyst 
layer  (i.e.,  as  a  region  or  interface)  was  not  significant.  For  de¬ 
tails  of  the  transport  phenomena,  the  treatment  of  the  catalyst 
layer  would  be  more  significant.  The  single-domain  simply 
used  a  finer  computational  grid  around  the  catalyst  layers. 
This  may  have  accounted  for  the  multi-domain  solution  re¬ 
quiring  extra  time  since  the  Butler- Volmer  equation  entails 
more  computationally  expensive  calculations.  So,  although 
in  this  case  the  multi-domain  solution  required  more  compu¬ 
tational  time,  it  gave  more  precise  information  about  species 
concentration  in  the  catalyst  layers.  The  actual  convergence 
time  for  any  model  depends  on  the  degree  of  accuracy  and  the 
efficiency  of  the  solution  algorithm.  It  cannot  be  concluded 
that  the  single-domain  method  is  faster  than  the  multi-domain 
method. 

For  the  one-dimensional  case,  the  solution  process  is  much 
simplified.  In  the  solutions  presented  in  this  paper,  the  cell 
current  density  was  given  and  the  cell  voltage  was  calcu¬ 
lated  for  that  particular  current  density.  For  two-  and  three- 
dimensional  cases,  it  is  difficult  to  specify  the  cell  current 
density  since  it  is  defined  as  the  area  integrated  average  across 
the  catalyst  layers.  In  such  cases,  the  cell  current  density  must 
be  specified,  but  an  iterative  approach  must  be  used  until  the 
area  integrated  average  current  density  across  the  catalyst 
layers  match  the  specified  cell  current  density.  Then  the  local 
overpotentials  can  be  computed. 

The  most  significant  measure  of  performance  of  a  fuel 
cell  is  the  polarization  curve.  For  PEMFCs,  water  manage¬ 
ment  is  also  a  critical  issue.  Polarization  effects  can  be  stud¬ 
ied  accurately  using  one-dimensional  models.  Comprehen¬ 
sive  water  management  requires  2-D  or  3-D  modeling.  How¬ 
ever,  a  one-dimensional  model  can  approximate  the  required 
water  addition/removal  required  for  optimum  operation  of 
the  Nafion®  membrane,  i.e.,  maintaining  the  membrane  in 
a  fully  humidified  condition.  Results  presented  here  and  in 
the  literature  for  one-dimensional  models  show  good  agree¬ 
ment  between  theoretical  and  experimental  polarization  ef¬ 
fects.  However,  one-dimensional  models  cannot  accurately 
predict  the  spatial  concentration  and  fluxes  of  species  within 
the  fuel  cell.  Such  information  is  needed  for  the  effective 
design  and  loading  of  catalysts.  Two-dimensional  models 
can  accurately  predict  spatial  variations  in  species  concen¬ 
tration  and  fluxes  for  simple  flow  regimes,  e.g.,  straight  flow 
channels  (co-flow  and  counter-flow  arrangements).  How¬ 


ever,  for  more  complex  flow  regimes,  such  as  interdigi- 
tated  flow  and  serpentine  flow  channels,  three-dimensional 
modeling  is  required  to  accurately  describe  the  transport 
phenomena. 

7.  Conclusions 

A  review  of  recent  literature  in  PEM  fuel  cell  modeling 
was  presented.  Fuel  cell  models  can  be  categorized  as  ana¬ 
lytical,  semi-empirical  or  mechanistic.  Mechanistic  models 
can  be  further  subcategorized  based  on  the  solution  strategy, 
single-domain  or  multi-domain.  The  merits  and  demerits  of 
each  were  discussed. 

The  single-domain  is  less  cumbersome  in  that  no  inter¬ 
nal  boundary  conditions  and  conditions  of  continuity  need  to 
be  specified.  It  is  also  easier  to  incorporate  into  commercial 
CFD  codes.  As  a  result,  the  time  for  model  development  is 
shortened. 

The  model  of  Bernardi  and  Verbrugge  [11]  was  taken  as 
a  case  study  and  solved  using  both  approaches.  Both  models 
accurately  predicted  the  polarization  effects  and  water  man¬ 
agement  requirements. 
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